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Summary 

The  ray  solution  for  stationary  hydrostatic  mountain  waves  has  a  singularity  along  the  vertical  axis  directly 
over  the  mountain.  We  use  Maslov’ s  method  to  improve  the  ray  prediction.  The  ray  solution  is  determined  in  the 
wave-number  domain  and  is  then  mapped  by  inverse  Fourier  transform  to  give  a  spadal  description  of  the  wave 
field  that  approximates  the  linear  solution  over  the  mountain  and  elsewhere.  We  develop  Maslov’s  method  for  a 
medium  with  a  height-dependent  mean  wind  and  mean  buoyancy,  and  we  compare  it  with  the  traditional  transform 
method  and  with  the  results  of  a  numerical  simulation. 
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1.  Introduction 

In  a  previous  study,  we  developed  a  ray-tracing  scheme  for  mountain  waves  that 
includes  the  prediction  of  wave  amplitudes  along  the  ray  (Broutman  el  al.  2001).  We 
initialized  the  rays  with  a  near-field  stationary-phase  approximation,  and  we  used  the 
ray-tracing  to  examine  the  far-field  approach  to  critical  layers  in  three  dimensions  with 
rotation.  We  were  unable,  however,  to  analyse  the  region  directly  over  the  mountain. 
Ray  theory  is  inaccurate  there  because  of  the  breakdown  of  the  slowly  varying  approxi¬ 
mation. 

In  the  present  study  we  first  show  that  this  breakdown  is  associated  with  a  special 
type  of  caustic  that  cannot  be  treated  by  the  usual  correction  techniques,  such  as  an  Airy- 
function  matching.  We  then  use  Maslov’s  method  to  improve  the  ray  prediction.  The 
method  was  developed  in  the  1960's  by  V.  R  Maslov  as  a  short-wavelength  asymptotic 
technique  for  problems  in  quantum  mechanics.  An  introduction  to  the  method  has 
been  given  by  Maslov  and  Fedoriuk  (1981),  Kravtsov  and  Orlov  (1993),  Thomson  and 
Chapman  (1985),  and  Ziolkowski  and  Deschamps  (1984).  Another  helpful  reference,  on 
which  our  own  work  is  most  closely  based,  is  Brown  (2000). 

To  derive  Maslov’s  solution,  we  express  the  ray  solution  as  a  function  of  wave- 
number  coordinates,  rather  than  in  the  usual  way  as  a  function  of  spatial  coordinates. 
In  wave-number  coordinates,  the  rays  in  our  model  remain  well  separated  and  do  not 
form  caustics.  The  ray  solution  is  then  mapped  from  the  wave-number  coordinates  to 
the  spatial  coordinates  by  an  inverse  Fourier  transform.  The  result  reproduces  the  usual 
spatial  ray  solution  where  that  ray  solution  is  accurate,  but  also  corrects  the  caustic 
singularity  in  the  region  over  the  mountain. 

There  is  another  way  to  derive  Maslov’s  solution,  by  a  variation  of  the  standard 
transform  technique.  The  difficulty  with  the  standard  technique  is  that  the  normal  modes 
are  required,  and  these  can  rarely  be  expressed  in  a  useful  analytical  form.  Maslov’s 
method  amounts  to  a  simplification  in  which  the  normal  modes  are  approximated  by  ray 
theory.  Shutts  (1998)  has  made  this  simplification  in  an  analysis  of  hydrostatic  mountain 
waves  in  a  turning  wind.  We  show  in  section  3(c)  that  Maslov’s  method  gives  Shutts’s 
solution  as  a  special  case.  The  consideration  of  normal  modes  also  leads  to  a  criterion 
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for  the  validity  of  Maslov’s  solution,  which  depends  on  the  largeness  of  the  Richardson 
number.  Note  that  the  standard  ray  method  fails  over  the  mountain  for  all  Richardson 
numbers. 

In  section  4,  we  compare  Maslov’s  method  with  a  numerical  simulation  in  Shutts 
and  Gadian  (1999),  finding  reasonable  agreement.  They  study  mountain  waves  in  a 
turning  wind  using  a  nonlinear-  non-hydrostatic  model,  though  their  choice  of  a  small 
wide  hill  for  this  simulation  favours  a  linear  hydrostatic  response. 

Throughout  this  paper  we  assume  linear  stationary  hydrostatic  waves  without 
rotation.  The  medium  is  horizontally  homogeneous,  with  smooth  height-dependent 
mean  wind  and  mean  buoyancy.  These  assumptions  simplify  the  implementation  of 
Maslov’s  method  in  significant  ways.  For  example,  the  combination  of  horizontal 
homogeneity  and  hydrostatic  waves  prevents  caustics  from  occurring  in  the  wave- 
number  coordinates.  Maslov’s  method  handles  cases  in  which  there  are  caustics  in  both 
spatial  and  wave-number  coordinates,  but  the  implementation  is  more  difficult  (e.g. 
Brown  2000).  Other  potential  difficulties  associated  with  the  inclusion  of  rotation  and 
non-hydrostatic  effects  are  noted  at  the  end  of  section  3. 

We  use  the  Boussinesq  approximation  until  section  4.  In  that  section,  Maslov’s 
solution  is  scaled  in  the  standard  way  to  account  for  the  decrease  in  mean  density  with 
height  (Gill  1982,  section  6.14).  This  is  necessary  for  the  comparison  with  the  simulation 
by  Shutts  and  Gadian  (1999),  whose  model  uses  the  anelastic  form  of  the  continuity 
equation. 


2.  Ray  theory  for  stationary  hydrostatic  mountain  waves 

We  first  describe  the  traditional  ray  solution  for  three-dimensional  stationary  hydro¬ 
static  mountain  waves.  Examples  are  given  by  Smith  (1980)  for  the  case  of  a  uniform 
medium,  and  by  Shutts  (1998)  and  Broad  (1999)  for  cases  with  wind  shear.  Our  aim  here 
is  to  understand  why  ray  theory  fails  directly  over  the  mountain.  We  use  coordinates 
x  —  (x,  y)  in  the  horizontal  directions  and  z  in  the  vertical  direction.  The  mean  wind 
is  horizontal  and  has  the  form  U  =  (U (z),  V (z)),  and  the  mean  buoyancy  frequency  is 

N(z). 

The  mountain  waves  have  vertical  wave  number  m  and  horizontal  wave-number 
vector 

k  =  (k,  l)  =  K(cos  (f>,  sin  <fi),  (1) 

with  magnitude  K  and  angle  <fi  from  the  positive  x-axis.  We  take  the  positive  x-axis  to 
point  in  the  downwind  direction  at  ground  level.  Then  the  combination  k  <  0,  m  <  0 
and  intrinsic  frequency  co  >  0  gives  stationary  waves  with  upward  group  propagation. 
The  intrinsic  frequency  and  wave  number  satisfy 

co  =  -KN/m  =  -kU  -IV.  (2) 

The  first  relation  in  (2)  is  the  hydrostatic  form  of  the  internal-wave  dispersion  relation 
under  present  assumptions,  while  the  second  relation  is  the  condition  for  stationary 
waves  which,  by  definition,  have  zero  frequency  relative  to  the  ground.  The  group 
velocity  relative  to  the  ground  has  horizontal  components 

cg  =  (cgU  cgi)  —  (U  —  kN /Km,  V  —  IN /Km)  (3) 

and  vertical  component  cg3  =  KN /mr . 

It  is  convenient  to  use  z  as  the  independent  ray  variable,  following  Shutts  (1998) 
and  Broad  (1999).  The  ray  path  x(z)  is  then  determined  from 

dx/dz  —  Cg/cg3. 


(4) 
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All  rays  are  launched  from  the  origin,  consistent  with  a  stationary-phase  analysis 
(Broutman  et  al.  2001).  Since  the  medium  is  horizontally  homogeneous,  k  and  1  are 
constant  along  the  ray. 

It  follows  from  (3)  and  (2)  that 


k  •  cg  =  0.  (5) 

It  follows  from  (4),  with  the  initial  ray  position  at  the  origin,  from  (5),  and  from  the 
constancy  of  k  along  the  ray,  that 

k#x  =  0.  (6) 

Therefore 

x  ||  cg.  (7) 

This  shows  that  the  horizontal  group  velocity  vector  always  lies  on  a  line  extending 
radially  from  the  origin.  Each  ray  is  thus  constrained  to  remain  in  a  single  vertical 
plane  throughout  its  propagation.  The  propagation  plane,  with  equation  x/y  =  —l/k  — 
—tan  <f>,  passes  through  the  origin  and  is  normal  to  (k,  TO).  All  rays  with  the  same  value 
of  cp  share  the  same  propagation  plane.  Similar  points  have  been  made  by  Shutts  (1998) 
and  by  Broad  (1999).  Note  that  the  conditions  (5 ) — (7 )  hold  only  in  the  hydrostatic  limit 
without  rotation. 

We  illustrate  this  behaviour  with  an  example  taken  from  the  model  of  Shutts  (1998). 
He  considers  the  wind  profile 

U(z)  =  (t/0,  Az),  (8) 

where  Uq  and  A  are  constants  and  z  —  0  is  the  ground.  Shutts’s  model  reduces  to  that 
of  Smith  (1980)  when  A  =  0.  In  the  example  plotted  below,  we  set  A  =  0.003  s-1, 
Uq  —  10  m  s_1,  and  N  —  0.01  s-1,  values  used  in  Shutts’s  calculations. 

Two  sets  of  five  ray  paths  each  are  shown  in  Figs.  1(a)  and  (b).  The  mean- 
wind  profile  (8)  used  in  these  calculations  is  indicated  in  Fig.  1(a).  Also  indicated  by 
the  dashed  line  in  each  plot  is  the  projection  of  the  ray  paths  onto  the  z  =  0  plane. 

The  projection  is  a  single  straight  line  because  the  rays  in  each  set  have  the  same  value 

of  (j)  and,  therefore,  share  the  same  propagation  plane.  In  terms  of  the  polar  coordinate 
angle  9  —  tan-  1  (y/x),  the  propagation  planes  in  Figs.  1(a)  and  (b)  are  defined  by 
6  =  50°  and  9  —  —50°  respectively.  In  our  wave-number  notation  (1)  this  corresponds 
to  wave-number  angles  <p  =  140°  and  <p  =  —140°,  respectively.  The  five  rays  in  each 
plot  have  K  values  equally  spaced  between  0.2/ L  and  1/L,  where  L  —  20  km  is  the 
mountain  width  used  by  Shutts  (1998). 

Rays  that  propagate  over  the  first  quadrant  (Fig.  1(a))  are  refracted  toward  decreas¬ 
ing  intrinsic  frequency.  Each  ray  asymptotes  toward  its  own  critical-layer  height.  This 
case  was  treated  in  detail  by  Shutts  (1998)  without  rotation  and  by  Broutman  et  al. 
(2001)  with  rotation.  These  rays  never  intersect  each  other  after  leaving  the  origin. 

Rays  that  initially  propagate  over  the  fourth  quadrant  (Fig.  1(b))  are  refracted 
toward  increasing  intrinsic  frequency.  Their  vertical  group  velocity  also  increases  with 
height,  and  they  encounter  an  increasingly  strong  cross-wind  V (z).  Eventually,  these 
rays  reverse  their  horizontal  direction  of  propagation,  moving  upwind  in  x  and  down¬ 
wind  in  y.  The  rays  loop  back  over  the  origin  and  intersect  at  a  single  point  on  the  vertical 
axis.  Rays  constrained  to  other  planes  with  —jz/2  <9  <  0  focus  at  other  heights  on  the 
vertical  axis,  forming  a  line  of  ray  singularities.  This  is  why  Shutts’s  stationary-phase 
solution  diverges  over  the  mountain  (see  Eq.  (67)  of  Shutts  1998). 

We  can  prove  the  general  result  illustrated  in  Fig.  1(b),  that  all  rays  in  the  same 
propagation  plane  intersect  the  vertical  axis  at  the  same  point.  We  return  to  the  case  of 
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(a)  (b) 


Figure  1.  Ray  paths  for  the  model  of  Shutts  (1998).  The  axes  are  labelled  in  km.  with  the  mean  wind  profile 
indicated  in  panel  (a)  only.  Rays  are  terminated  after  80  buoyancy  periods  of  propagation.  The  projection  of  the 
ray  paths  onto  the  x-y  plane  is  represented  by  the  dashed  line.  Two  types  of  behaviour  are  illustrated,  depending 
on  whether  the  rays  refract  toward  lower  intrinsic  frequency  or  toward  higher  intrinsic  frequency.  In  the  former 
case,  panel  (a),  each  ray  approaches  its  own  critical  level.  The  propagation  plane  is  0  =  tan~'(v/.r)  =  50°.  In  the 
latter  case,  panel  (b),  the  rays  loop  over  the  mountain  and  intersect  at  one  point  on  the  ;-axis.  The  propagation 

plane  is  d  =  —50°. 


arbitrary  U (z),  V (z)  and  N  (z).  It  is  sufficient  here  to  work  only  with  the  x -coordinate 
of  the  ray.  When  the  x  -coordinate  vanishes,  so  does  the  y-coordinate,  and  the  ray  has 
intersected  the  vertical  axis.  The  ray  equation  (4)  for  x(z)  becomes 


1=  [  ^-dz 
j  cg3 

_  tan  <t>  f  U (z)  tan  cl)-V(z) 

K  J  (C/(z)  +  V(z)  tan  0)2  ' 


(9) 


The  integral  in  (9)  is  taken  along  the  ray,  and  the  ray  is  launched  from  the  origin.  The  ray 
position  v  depends  in  general  on  K  but,  on  the  vertical  axis,  *  =  0  and  we  are  left  with  a 
relation  for  z(4> )  in  which  the  K  factor  cancels  out.  Thus  all  rays  in  the  same  propagation 
plane  (determined  by  0)  intersect  each  other  on  the  vertical  axis  at  precisely  the  same 
point  (or  set  of  points). 

The  curve  or  surface  on  which  neighbouring  rays  intersect  each  other  is  called  a 
caustic  (Lighthill  1978;  Kravtsov  and  Orlov  1993;  Sonmor  and  Klaassen  2000).  Each 
neighbouring  ray  has  a  slightly  different  wave  number,  so  the  wave-number  gradient 
diverges  on  the  caustic  in  violation  of  the  slowly  varying  approximation.  The  ray 
prediction  for  the  wave-action  density  also  diverges,  depending  inversely  on  the  size 
of  a  volume  element  enclosed  by  the  neighbouring  rays. 

There  are  standard  ways  to  correct  the  ray  solution  near  certain  types  of  caustics,  but 
this  caustic  is  not  one  of  those  types.  The  simplest  and  most  familiar  caustic  is  corrected 
with  an  Airy  function  (Lighthill  1978)  and  is  characterized  by  the  intersection  of  two 
neighbouring  rays  at  each  point  on  the  caustic.  More  complicated  caustics  require  differ¬ 
ent  correction  functions  and  involve  the  intersection  of  more  rays.  Three  neighbouring 
rays  intersect  at  the  tip  of  a  cusped  caustic,  for  example,  whose  correction  function  can 
be  found  in  the  monograph  by  Kravtsov  and  Orlov  (1993). 

As  we  have  just  shown,  the  caustic  over  the  mountain  is  intersected  by  an  infinite 
number  of  rays  at  each  point.  Such  caustics  arise  only  in  very  special  circumstances  and 
are  structurally  unstable  (Berry  198 l,p.  480).  That  is,  they  generally  break  apart  with  a 
small  perturbation  to  the  problem.  A  minor  non-hydrostatic  correction,  for  instance,  is 
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enough  to  eliminate  the  mountain-wave  caustic  entirely.  By  contrast,  the  simple  Airy- 
function  caustic  is  one  of  several  kinds  of  caustics  (including  the  cusped  caustic)  that 
survive  perturbation  and  are  structurally  stable.  A  slight  perturbation  to  a  problem  with 
a  stable  caustic  merely  leads,  in  general,  to  a  slight  movement  in  the  caustic’s  position. 

Broutman  et  al.  (2001)  computed  non-hydrostatic  ray  solutions  for  the  Shutts 
(1998)  model.  Although  the  non-hydrostatic  ray  paths  do  not  quite  intersect  over  the 
mountain,  they  still  focus  strongly  there,  causing  a  breakdown  of  the  slowly  varying 
assumption  and  erroneously  large  ray  predictions  for  the  vertical  displacement.  Another 
way  is  needed  to  correct  the  ray  solution,  and  for  this  purpose  we  now  turn  to  Maslov’s 
method. 


3.  Maslov’s  method 

Maslov’s  solution  is  the  inverse  Fourier  transform  of  the  ray  solution  expressed  in 
wave-number  coordinates.  To  derive  Maslov’s  solution,  we  start  with  the  more  familiar 
ray  solution  in  spatial  coordinates,  map  it  to  wave-number  coordinates  using  the  local 
ray  relation  between  wave  number  and  position,  apply  boundary  conditions,  and  then 
map  back  to  the  spatial  coordinates  by  inverse  Fourier  transform. 

For  stationary  mountain  waves  there  is  no  time  dependence,  and  the  z -dependence 
can  be  treated  parametrically.  A  ray  variable  describing  the  waves,  such  as  the  vertical 
displacement,  is  written  in  spatial  coordinates  as 

?;(x,  z)  =  a(x,  z )  e1^*’*)  (10) 

and  in  wave-number  coordinates  as 

rj(k,  z)  —  a(k,  z)  e1^’^.  (11) 

As  before  x  =  (x ,  y)  and  k  =  (k,  l).  The  inverse  Fourier  transform  of  rj( k,  z)  gives 
Maslov’s  solution 

/OO  PO O 

/  rj( k,  z)  elk‘x  d k  d /.  (12) 

-OO  J  — OO 

(a)  The  spatial  domain 

The  vertical-displacement  amplitude  |a(x,  z)  \  in  (10)  can  be  obtained  in  ray  theory 
from  conservation  of  wave-action.  The  wave-action  density  A  (x,  z)  and  the  wave-energy 
density  E  (x,  z)  are  related  by  A  =  E fed,  where  cd  is  the  intrinsic  frequency.  For  gravity 
waves  without  rotation 

A  =  ±p0\a\2N2/ad,  (13) 

where  po  is  the  mean  density.  In  a  steady-state  application  such  as  ours,  A  satisfies 

V  •  (cgA)  +  3(cg3A)/9z  =  0,  (14) 

where  V  =  (d/dx,  d/d y).  Integration  of  (14)  over  a  ray-tube  volume,  followed  by  an 
application  of  the  divergence  theorem  leads  to  (cf.  Broad  1999) 

/  x,  y  \ 

cg3  A  J  I  - —  I  =  constant  (15) 

W  yoJ 

following  the  ray.  Here  xo,  >'o  is  an  initial  point  on  the  ray  corresponding  to  z  =  zo- 
The  Jacobian  in  (15),  defined  by 

x,  y  \  3x  dy  dx  dy 

x0,  yoJ  dxo  dyo  dy0dx0’ 


J 


(16) 
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is  evaluated  at  constant  z  and  accounts  for  the  change  in  A  due  to  the  horizontal 
spreading  of  rays.  The  factor  of  cg 3  in  (15)  accounts  for  the  change  in  A  due  to  the 
vertical  spreading  of  rays. 

The  phase  1//  in  the  spatial  solution  (10)  is  given  by 


VAx,  z) 


k(x,  z) '  x  + 


m(x,  z! )  dz'. 


(17) 


All  rays  originate  from  the  origin,  where  t//  =  0,  consistent  with  the  stationary-phase 
solution  (Shutts  1998;  Broutman  et  al.  2001).  For  the  hydrostatic  case,  we  noted  in 
section  2  that  k  •  x  vanishes.  We  retain  this  term  for  now.  It  will  cancel  out  below  (in 
(24))  without  recourse  to  the  hydrostatic  approximation. 


( b )  The  wave-number  domain 

Next  we  determine  the  wave-number  solution  7f(k,  z)  in  (1 1).  Stalling  with  conser¬ 
vation  of  wave-action  (15),  we  re-express  the  Jacobian  using  the  chain  rule  to  yield 

c«3-4,(rr),(i77;)=cons“t  <18) 

The  product  of  A  with  the  first  Jacobian  in  ( 1 8)  transforms  the  wave-action  from  a  spatial 
density  to  a  wave-number  density  A(k,  z),  i.e., 


A(k,  z) 


A(x,  z)  J 


(19) 


The  second  Jacobian  in  (1 8)  is  constant  along  the  ray  in  our  problem.  Thus 

cg3A  =  constant.  (20) 

In  the  wave-number  domain,  the  vertical  flux  of  wave-action  is  constant  following 
the  ray.  Figure  2  helps  to  illustrate  this  point.  The  left  panel  shows  rays  in  the  spatial 
domain.  They  diverge  strongly  in  all  three  dimensions  as  they  move  away  from  the 
mountain  centred  at  the  origin.  The  right  panel  shows  the  corresponding  rays  in  the 
wave-number  domain.  Being  parallel  to  the  .“-direction,  they  diverge  or  converge  only 
with  respect  to  the  2-direction.  This  implies  the  form  of  wave-action  conservation  given 
in  (20). 

Ray  initialization  is  simple  in  the  wave-number  domain.  The  rays  can  be  initialized 
directly  from  the  mountain  profile  rather  than  through  the  near- Held  stationary-phase 
analysis  used  by  Broutman  et  al.  (2001)  to  initialize  rays  in  the  spatial  domain.  The  ver¬ 
tical  displacement  rj( k,  z)  in  (11)  is  equated  at  z  =  0  to  the  Fourier  transform  h( k)  of 
the  mountain  profile  h(x).  Since  1 Jr  vanishes  at  z  =  0  (see  (25)  below),  we  are  left  with 

a(k,  0)  =  h(k).  (21) 

Next  we  relate  A  to  rj  in  (1 1)  by 

A  =  \po\d\2N2 /d).  (22) 

analogous  to  the  definition  in  the  spatial  domain  (13)  and  analogous  to  the  plane- 
wave  polarization  relations  between  Fourier  amplitudes  of  different  variables  (Gill  1982, 
section  6.5).  Combining  (20),  (21),  and  (22)  yields 


a( k,  z) 


h(  k) 


~cg3(k,  0)  S(k,  z)l1/2  N(Q) 
_cg3(k,  z)  m(k,  0).  N(z) 


(23) 
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(a)  (b) 


Figure  2.  Rays  near  the  mountain  in  (a)  the  spatial  domain  and  (b)  the  wave-number  domain.  The  length  of  the 
ray  vectors  is  not  significant.  The  mountain  is  centred  at  the  origin  in  (a).  The  spatial  coordinates  are  labelled 
in  km.  The  wave-number  coordinates  k,  l  are  normalized  by  the  mountain  width  L  =  20  km.  We  have  ignored 
refraction  because  the  dominant  effect  near  the  mountain  is  geometrical  spreading.  Refraction  by  vertical  shear 
causes  the  rays  to  curve  in  the  spatial  domain  but  not  in  the  wave-number  domain. 


The  phase  in  the  wave-number  domain  x[r  is  related  to  the  phase  in  the  spatial- 
domain  x[r  by  the  Legendre  transformation  (e.g.  Brown  2000;  Ziolkowski  and 
Deschamps  1984) 

V/(k,  z)  =  VKx,  z)  —  k  •  x(k,  z).  (24) 

Substituting  from  (17)  gives  a  cancellation  of  k  •  x  terms,  leaving 

V^(k,  z)—f  m(k,  z!)  dz'.  (25) 

J  o 

The  cancellation  of  the  k  •  x  terms  occurs  only  because  k  is  constant  along  the  ray.  More 
generally  t//  depends  on  /  k  •  dx,  whereas  the  term  in  the  Legendre  transformation  is 
always  k  •  x. 


(c)  The  Maslov  solution 

Substituting  the  amplitude  solution  (23)  and  the  phase  solution  (25)  into  the  Fourier 
integral  (12)  gives  the  Maslov  solution 


pm(x,  z)  = 


ATO) 

N(z) 


OC 

-( 


H  k) 


cg3(k,  0)  S(k,  z) 
Lcg3(k,  z)  S(k,  0)  J 


1/2 


;i/o"«(k.-')dT  eik*x  d/t  d/. 


(26) 

In  the  case  of  a  uniform  medium,  the  ratio  in  square  brackets  is  unity,  as  is  the  ratio 
of  N  values,  and  /()  m(k,  z!)  dz'  =  m(k)z.  Thus  (26)  reduces  to  the  integral  solution  in 
the  uniform-medium  model  of  Smith  (1980,  Eq.  (11)).  For  constant  N  and  hydrostatic 
waves,  (26)  simplifies  to 


pm(x,  z)  = 


/oo  no 

-oo  J  — n 


H  k) 


m(k,  z) 

L  m  (k,  0)  J 


1/2 


*  Jo  dz'  eik-x  ^  d/_ 


(27) 


We  now  show  that  (27)  agrees  with  the  solution  derived  by  Shutts  (1998)  for  a 
particular  mean- wind  profile.  This  serves  as  a  check  on  our  derivation  and  illustrates  the 
relation  between  Maslov’s  method  and  a  standard  transform  method  involving  normal 
modes.  It  also  leads  to  a  condition  for  assessing  the  validity  of  Maslov’s  solution. 
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Shutts  (1998)  examines  stationary  hydrostatic  mountain  waves  in  a  medium  with 
constant  buoyancy  frequency  N  and  with  mean  wind  velocity  (U o,  A z),  where  Uq  and 
A  are  positive  constants  and  z  =  ()  is  the  ground.  Shutts  derives  an  equation  for  the 
normal  modes  in  terms  of  the  coordinate  £  defined  by 

1 

£  = - ,  (28) 

1  +  (Az  tan  4>)/ Uq 


where  4>  =  tan  1  l/k.  (Shutts’s  results  are  given  here  in  his  notation,  which  differs  from 
ours  in  that  his  k  and  m  are  positive.)  The  normal-mode  equation  is 


3  2rf 

H2 


/x  ~ 

+  p-i?  =  o, 


(29) 


where  /x  =  N/(A  sin  (f>).  As  noted  by  Shutts,  the  exact  solution  to  (29)  has  the  form  (  y 
with  v  —  1/2  it  i(/x2  —  1/4)1/2. 

At  this  point  in  his  derivation,  Shutts  makes  two  assumptions  that  amount  to 
invoking  the  ray  approximation.  First,  he  assumes  that  the  Richardson  number  N  2/ A2  is 
large,  so  that  /x2  3>  1/4.  Secondly,  he  considers  only  waves  with  upward  group  velocity, 
corresponding  to  the  negative  root  in  the  above  expression  for  v.  This  ignores  partial 
reflection,  which  is  also  ignored  in  standard  ray  theory.  After  applying  the  boundary 
condition  at  z  —  0,  Shutts  finds  that 

?j  =  £1/2A(  k)e_i/ilogf.  (30) 


This  result  can  also  be  derived  from  (29)  by  using  a  standard  formula  for  the  WKBJ 
approximation  (e.g.  Gill  1980,  Eq.  (8.12.7)). 

Shutts’s  solution  in  the  spatial  domain  is  the  two-dimensional  Fourier  integral  of 
(30).  But  this  is  exactly  Maslov’s  solution  (27).  To  show  that  the  two  are  equivalent, 
we  first  note  from  (28)  that  for  each  normal  mode  £  =  <7>(k,  0)/m(k,  z),  so  that  in 
the  hydrostatic  approximation  £  =m(k,  z)/m(k,  0).  It  can  also  be  shown  that  the 
z-derivative  of  the  phase  function  — /x  log  £  in  (30)  is  equal  to  m.  Thus  (30)  can  be 
rewritten  as 


rj  =  h  (k) 


m(k,  z ) 

|_m(k,  0)  J 


1/2 


d  Jo  m(kz')dz' 


(31) 


which  is  the  quantity  in  Maslov’s  Fourier  integral  (27). 

Maslov’s  method  is  thus  a  simplification  of  the  standard  transform  method  involv¬ 
ing  normal  modes.  Instead  of  solving  for  the  normal  modes  exactly,  Maslov’s  method 
uses  their  ray  approximation.  The  exact  normal  modes  cannot  generally  be  determined 
in  a  useful  analytical  form,  but  the  ray  approximation  is  a  simple  expression  that  is  easily 
calculated. 

The  accuracy  of  Maslov’s  solution  (27)  thus  depends  on  how  well  the  ray  approxi¬ 
mation  represents  the  normal  modes.  The  normal  modes  are  z -dependent,  so  we  expect 
the  ray  approximation  to  be  valid  when  the  vertical  wave  number  is  slowly  varying,  i.e. 
when  |m_23m/3z|  1  (see  Lighthill  1978,  p.  324). 

We  rewrite  the  dispersion  relation  (2)  as 

KN  N 

m  — - = - .  (32) 

kU  +  IV  U  cos  4>  +  V  sin  (p 

We  allow  U  and  V  to  be  general  functions  of  z,  but  keep  N  constant  (following  Shutts 
(1998)  and  as  in  our  comparison  below  with  Shutts  and  Gadian  (1999)).  We  then  obtain 

1  dm  Uz  cos  <p  +  Vz  sin  <p 
m2  3z 


N 


(33) 
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Hence  the  ray  approximation  to  the  normal  modes,  and  consequently  Maslov’s  method, 
is  accurate  for  our  model  when  the  Richardson  number  N2 /(U:  +  V2)  is  large. 

The  condition  |m~23m/3z|  1  used  here  should  not  be  confused  with  the  con¬ 

dition  for  ray-theory  validity  in  the  spatial  domain.  We  do  not  know  the  appropriate 
condition  in  the  spatial  domain,  but  it  cannot  be  \m~2dm/dz\  1.  This  is  obvious  from 
the  discussion  in  section  2  and  from  the  stationary-phase  calculation  of  Shutts  (1998); 
ray  theory  breaks  down  in  the  region  over  the  mountain  even  if  the  Richardson  number 
is  large  (or  infinite,  as  in  the  case  of  Smith  1980).  Over  the  mountain,  there  are  strong 
gradients  in  the  horizontal  wave  numbers,  which  violate  the  slowly  varying  assumption. 
These  strong  gradients  are  irrelevant  when  seeking  a  ray  approximation  to  the  normal¬ 
mode  equation  because  the  horizontal  wave  numbers  are  then  heated  as  constants. 

(Note  that  in  the  derivation  of  (33)  from  (32),  <p  was  not  differentiated  with  respect 
to  z  because  k  and  /,  and  hence  <fi,  are  fixed  for  each  mode.  For  the  spatial  domain, 
k  and  1  are  considered  to  be  functions  of  all  spatial  coordinates.  But,  as  discussed 
in  section  2,  each  vertical  plane  contains  waves  of  the  same  4>,  so  <pz  =  0.  Thus  the 
expression  given  in  (33)  is  correct  for  the  spatial  domain  as  well.  It  is  the  validity 
condition  |m_23m/3z|  1  that  is  not  appropriate  for  the  spatial  domain.) 

In  deriving  Maslov’s  solution  (26),  we  did  not  make  the  hydrostatic  approximation. 
There  is  a  complication,  however,  in  applying  Maslov’s  theory  to  non-hydrostatic  waves. 
The  ray  approximation  to  the  normal  mode  diverges  at  a  buoyancy-frequency  turning 
point  where  cg 3  vanishes.  The  divergence  is  evident  from  the  constancy  of  cg3  A(k,  z)  in 
(20)  and  represents  a  caustic  in  the  wave-number  domain.  One  possible  remedy  is  to  use 
an  Aiiy -function  approximation  to  the  normal  mode  near  its  caustic,  though  we  have  not 
tried  this,  nor  have  we  yet  considered  trapped  modes. 

We  can  add  rotation  to  the  Maslov  solution,  but  again  there  are  complications.  In  the 
rotating  case,  another  factor  is  introduced  into  the  relation  between  vertical  displacement 
and  wave-action  density  in  (13)  (see  Eq.  (34)  of  Broutman  et  al.  (2001)).  This  factor 
carries  through  the  subsequent  analysis  and  appears  in  the  integrand  of  (26)  as  another 
ratio  between  values  at  z  =  0  and  at  z.  The  first  complication  (e.g.  Inverarity  and  Shutts 
2000,  p.  2712)  is  that  in  the  general  rotating  case  the  mean  density  has  horizontal 
variations  imposed  by  the  thermal- wind  relation.  Thus  the  horizontal  wave  numbers  k,  1 
are  not  constant  along  the  ray,  though  such  variations  in  k,  1  may  be  unimportant  in  some 
cases  of  interest.  A  second  possible  complication  is  that  the  standard  ray  approximation 
to  the  normal-mode  equation  fails  to  satisfy  the  slowly  varying  assumption  near  a 
rotating  critical  layer  (Yamanaka  and  Tanaka  1984).  A  third  complication  is  the  presence 
of  an  infinitely  long  lee-wave  train  of  near- inertial  oscillations  (Shutts  2001).  This  causes 
errors  due  to  periodic  wrap-around  effects  when  the  Fourier  integral  in  Maslov’s  solution 
is  approximated  by  a  discrete  transform  (as  in  the  upcoming  example).  The  use  of  a 
complex  frequency  to  damp  the  inertia  oscillations  and  limit  their  horizontal  extent  is 
discussed  in  Shutts  (2001). 

In  the  following  we  continue  to  restrict  attention  to  the  hydrostatic  case  without 
rotation. 


4.  An  EXAMPLE  from  Shutts  and  Gadian  ( 1  999) 

We  compare  Maslov’s  method  with  the  numerical  results  shown  in  Fig.  1 1  of  Shutts 
and  Gadian  (1999).  Although  they  use  a  nonlinear  non-hydrostatic  model,  their  choice  of 
parameters  favours  a  linear  hydrostatic  response.  For  this  case,  the  response  compares 
well  with  the  linear  hydrostatic  theory  also  developed  by  Shutts  and  Gadian.  Normal 
modes  are  used  in  the  theory  but,  since  these  have  a  very  complicated  analytical  form. 
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Shutts  and  Gadian  found  it  simpler  to  integrate  the  normal-mode  equation  numerically, 
in  its  exact  form.  As  we  have  discussed  in  the  previous  section,  Maslov’s  method  uses  a 
ray  approximation  to  the  normal  modes. 

Shutts  and  Gadian  (1999)  have,  in  fact,  supplied  a  test  of  Maslov’s  method.  One  of 
their  numerical  simulations  is  compared  against  the  theory  of  Shutts  (1998)  who,  as  we 
mentioned  in  section  3(c),  derived  a  special  case  of  Maslov’s  solution.  The  agreement 
between  the  numerics  and  Shutts’s  theory  is  very  good  in  this  case,  as  indicated  by 
Figs.  8(d)  and  9  of  Shutts  and  Gadian. 

For  the  case  in  Fig.  1 1  of  Shutts  and  Gadian  (1999),  the  wind  profile  is  given  by 

U(z)  =  U o fcos(7T z / H),  sin (nz/H)],  (34) 

where  Uq  =  10  m  s  "1  and  H  =  12  km.  The  mountain  height  is 

Hr)  =  h0/(l  +  r2/L2)3'2,  (35) 

where  r  =  (x2  +  y2)1/2,  L  =  20  km,  and  ho  —  100  m.  In  terms  of  the  horizontal-wave- 
number  magnitude  K,  the  Fourier  transform  of  h  (r)  is  (Smith  1980;  Shutts  1998) 

h(K)  —  hoL2  exp (-KL).  (36) 


From  (34),  the  term  under  the  square  root  in  the  integrand  of  Maslov’s  integral  (27) 
becomes 


m(k,  z)/m( k,  0)  =  k/[k  cos(nz/H )  +  /  sin(jrz///)].  (37) 


The  phase  function  in  (27)  can  be  evaluated  analytically  as 


f 

J  o 


m(k,  z')  dz' 


2 NH 

TTtT 


{ k  tan(jz z/2H)  —  l\ 
-  +  tanh 

V  K  ) 


(38) 


The  mean-density  profile  in  Shutts  and  Gadian  corresponds  to  N  =  0.0113  s 
and  the  Richardson  number  N2 / |UZ  |2  is  approximately  19.  To  account  for  the  growth  in 
wave  amplitude  due  to  the  decrease  in  mean  density  with  height,  the  dependent  variables 
calculated  below  are  multiplied  by  exp(z/2//p)  (see  section  6.14  of  Gill  (1982)).  We  use 
a  density  scale  height  of  Hp  =  7.5  km,  consistent  with  typical  values  for  the  troposphere 
(Gill  1982,  p.  50). 

The  Maslov  solutions  for  the  x  and  z  perturbation  velocities,  u  and  w,  respectively, 
are  needed  in  order  to  compare  with  Shutts  and  Gadian’s  results.  The  conversion  from 
vertical  displacement  i]  in  (27)  to  u  and  w  is  accomplished  by  ray  approximation. 
For  gravity  waves,  the  required  ray  relations  are  w  =  — i corf  and  u  =  — kmw/K 2 .  The 
Maslov  expression  for  w,  for  example,  is  (27)  with  the  integrand  multiplied  by  —  ico. 

Maslov’s  Fourier  integral  is  approximated  using  a  discrete  Fourier  transform, 
with  ±k  wave  numbers  included  as  complex  conjugate  solutions  (Smith  1980).  In  the 
example  below  we  used  a  256  by  256  grid  of  k,  l  wave  numbers  (including  complex 
conjugates),  with  maximum  wave-number  magnitudes  of  kmax  =  8 / L  and  /max  —  6/L, 
L  being  the  mountain  width.  Doubling  kmAX  and  Zmax  (with  the  same  total  of  25  62  wave 
numbers)  led  to  very  similar  results,  except  for  some  noticeable  periodic  wrap-around 
effects  in  the  case  of  Fig.  3(c).  Halving  kmAX  and  /max  led  to  significant  short-scale  errors, 
especially  in  the  cases  of  Figs.  3(c)  and  (d),  where  peak  values  decreased  by  about  30%. 

In  Fig.  3,  we  plot  the  Maslov  solutions  for  u  (panels  (a)  and  (c))  and  w  (panels  (b) 
and  (d))  that  correspond  to  Fig.  1 1  of  Shutts  and  Gadian.  These  are  given  at  z  =  2.1  km 
(panels  (a)  and  (b))  and  at  z  —  6.3  km  (panels  (c)  and  (d)).  For  the  k,  1  grid  we  have 
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Figure  3.  The  Maslov  solution  for  the  case  in  Fig.  11  of  Shutts  and  Gadian  (1999).  The  plots  show  the 
perturbation  velocities  u  and  w  at  the  heights  indicated  above  each  panel.  The  contour  interval  (Cl)  and  the 
minimum  and  maximum  values  of  the  variables  are  (a)  Cl  =  0.05,  (—0.38,  0.35),  (b)  Cl  =  0.004,  (—0.035, 
0.030),  (c)  Cl  =  0.05,  (-0.25,  0.26),  (d)  Cl  =  0.002  (-0.018,  0.012). 


chosen,  the  full  domain  is  about  2011  km  by  2680  km.  Beyond  these  distances  the 
solution  wraps  around  periodically. 

The  Maslov  solutions  are  in  reasonable  agreement  with  Shutts  and  Gadian  (1999). 
For  example,  for  the  w  solution  in  Fig.  3(b),  the  Shutts  and  Gadian  contours  show 
a  minimum  value  somewhere  between  —0.028  and  —0.032  and  a  maximum  value 
somewhere  between  0.028  and  0.032.  The  corresponding  minimum  and  maximum 
Maslov  values  are  approximately  —0.035  and  0.030.  The  ranges  for  the  other  Maslov 
solutions  are  given  in  the  caption. 


5.  Discussion 

When  the  slowly  varying  approximation  breaks  down  in  the  spatial  domain,  the  idea 
of  Maslov’s  method  is  to  look  instead  to  the  wave-number  domain  for  a  slowly  varying 
ray  solution.  A  spatial  solution  is  then  obtained  by  inverse  Fourier  transform,  giving 
Maslov’s  result  (27).  Away  from  the  caustic  over  the  mountain,  the  inverse  Fourier 
transform  takes  on  its  stationary -phase  asymptotic  limit  and  thus  duplicates  the  standard 
spatial  ray  solution  in  which  one  ray  makes  the  dominant  contribution  at  each  position. 
Near  the  caustic  over  the  mountain,  a  spectrum  of  rays  contributes  significantly  to  the 
wave  field  at  each  position,  and  this  too  is  automatically  represented  by  the  inverse 
Fourier  transform.  In  other  applications,  other  types  of  caustics  occur,  each  properly 
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represented  in  Maslov’s  method  by  the  inverse  Fourier  transform  in  the  appropriate 
asymptotic  limit. 

Several  features  of  our  model  have  made  the  implementation  Maslov’s  method 
very  straightforward.  At  the  end  of  section  3  we  noted  the  potential  difficulties  that 
we  avoided  by  our  assumption  of  hydrostatic  waves  (which  eliminates  turning  points 
and  trapped  modes)  and  from  our  neglect  of  rotation.  We  also  emphasize  that  many 
applications,  unlike  the  present  one,  will  lead  to  caustics  in  both  the  wave-number  and 
spatial  domains.  Maslov’s  method  must  then  combine  ray  solutions  from  both  domains. 
Brown  (2000)  demonstrates  how  this  can  be  done  in  a  model  of  surface  gravity  waves. 

The  present  work  represents  another  step  in  our  development  of  ray-tracing  schemes 
for  internal  gravity  waves.  A  previous  study  concerned  the  (spatial-domain)  initializa¬ 
tion  of  the  ray-tracing,  and  was  discussed  by  Broutman  et  al.  (2001).  This  initializa¬ 
tion  is  suitable  for  a  number  of  sufficiently  idealized  models  involving  radiation  from  a 
localized  source,  including  those  models  that  describe  internal  gravity  waves  generated 
by  flow  over  topography,  by  a  collapsing  mixed  region  (Biihler  et  al.  1999;  Biihler  and 
McIntyre  1999),  and  possibly  by  buoyancy  oscillations  in  clouds  (e.g.  Piani  et  al.  2000). 
An  important  need  now  is  a  practical  way  to  correct  the  spurious  amplification  associ¬ 
ated  with  the  approach  to  a  caustic,  or  the  near  miss  of  a  caustic.  Our  initial  experience 
with  Maslov’s  method  here  shows  that  it  can  be  helpful  in  this  regard.  Further  work, 
especially  with  non-hydrostatic  and  rotation  effects  included,  is  required  to  understand 
the  potential  of  Maslov’s  method  and  the  range  of  applications  for  which  it  can  bring 
useful  insights. 
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